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ABSTRACT 


A comparison of numerical methods utilized by the finite 
element technique for solving a nonlinear nuclear reactor 
dynamics problem was conducted. Using the Crank-Nicolson, 
DVOGER (Gear) and Implicit Gear methods, the results showed 
Pacuinplictesto bewehe Superior method investigated. This 
is based on the fact that all three methods yielded the same 
steady state solutions; but, the Implicit Gear method used 
Significantly less CPU time and comparable storage to Crank- 
Nicolson. This was particularly apparent as the degrees of 
freedom were increased. In addition, the transient solution 
Gueall cases was betcer than that obtained in Crank-Nicolson 
and compared favorably to that of Gear's method. 

The other noteworthy result was in the effect of the 
emEor Chiterion en solution. It was shown that for a range 


orton cue ty 


to 1.0, the steady state solution value 
remained the same. This results in a significant reduction 
in computer processing time since the time required decreases 


substantially as the error conditions imposed are relaxed. 
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Eo LNERODUCTION 


i PURPOSE 

This research project has been undertaken to compare 
several numerical methods of solving a nonlinear nuclear 
reactor dynamics problem. Three methods have been investi- 
gated in this thesis, including the Crank-Nicolson, the 
DVOGER (Gear) and the Implicit Gear methods of solution. 

A nuclear reactor dynamics problem with temperature- 
dependent feedback, when it entails either a non-homogeneous 
or multi-region reactor, results in a nonlinear field equa- 
tion in space and time. This problem does not lend itself 
to solution by analytical means [5, 6]. However, when the 
physical and neutronic properties of the problem are known, 
a model can be formulated using the finite element method 
which will yield the transient and steady state flux solu- 
tions. In particular, the partial differential equations 


investigated were of the form 


aot = bP? w+ cw -wh" (1) 


where a, b, ¢, and w are constants andw({r, z, t}) is the 
flux. The finite element method reduces Equation (1) to 


the system of ordinary differential equations 


N N 


Dare a 5 pee ic - 
ga Agy ¥j0t) = 52) Bay jt) i= 1, -.-,N (2) 


when the nonlinear term is linearized. Nguyen and Salinas [5] 
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and more recently Olsen [6] discuss the problem and methods 
ee solution. 

In this work, a comparison of three numerical integration 
schemes for Equation (2) was made. The comparisons will be 
made on several bases: computer storage requirenents, com- 
puter processing time, rapidity with which solution was 
obtained and the relative accuracy of the solution. 

At steady state, it is expected that all three methods 
Will provide the same solution value. This is due to the 
fact that at steady state the time derivative of the flux 
(v), Equation (2), is zero. 

Mponder toscxplore the relative value of the various 
methods, the model has been discretized into finite element 
grids of various degrees of freedom (DOF). The effect on 
the solution by finer discretizations of the finite element 
has been investigated to determine if the various methods of 
Polutioneare Similarly inflwenced to provide a better solution 
for a finer mesh. To have a method of solution relatively 
independent of mesh size would greatly reduce computer storage 
and processing time requirements since a larger grid with 
fewer elements could be utilized. 

In order to test the flexibility of the various equation 
solvers, an initial disturbance was introduced at different 
points in the model. This provided both a check of the 
ability of the method to accept a random disturbance as well 
as information on how rapidly a solution is obtained with a 


Pattteular disturbance input. Two nodal points of the system 


di 





feue considered, a point at the origin and a point on the 
core-reflector interface. This was done to determine if the 
tracking ability was consistent throughout the model. 

Additional comparisons investigated include the effects 
of the convergence criterion on the solution for all methods 
and the effects of the size of the time step on the Crank- 
Nicolson method. 

Modifications have been made to the programs provided in 
wee. (6) which was the basis of this research project. In 
Bit Wemk, COre properties were arbitrarily assigned to the 
reflector elements at the interface. This occurred because 
the material and nuclear properties were provided at the nodal 
points rather than by elements. The properties were intro- 
duced in this manner as a means of computer storage reduction. 
Regardless of the mesh size, there were always fewer nodes 
than elements. In this project, all properties were provided 
on an element basis to eliminate this discrepancy while at 
tne same time sacrificing the minor storage savings it repre- 


sented. 
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II. DATA GENERATORS 


meee GENERAL 

The data generators formulated in this work are useable 
only for regular rectangular grids (Figures 1-3). Having 
selected the number of horizontal and vertical points for 
the discretization of the model, the data generators will 
provide the numbering of each nodal point, the horizontal 
and vertical position of each node and the nodal neighbors 
of each node. 

In addition, the outer boundary nodes, at which the 
dymemre flux is zéro, will be numbered last in order to 
reduce the number of equations required by the particular 
method of solution being used. This will substantially 
Meercase thre Storage requirements. For example, in a one 
hundred thirty-two node discretization, only one hundred ten 


equations will be solved. 


B. PROPERTY INPUTS 
i EES eC 

A simple data generator has been provided which will 
produce a data deck containing the physical properties of 
piemredcton core and reflector in the format required by the 
Crank-Nicolson and DVOGER (Gear) methods. These properties 
are neutron velocity, neutron diffusion length, neutron 
absorption cross-section, neutron fission cross-section and 


reactivity temperature coefficient. 
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This generator will accommodate any size mesh and 
will process an indefinite number of grids simultaneously. 

The user must provide this routine with a data deck 
which contains the total number of elements in the grid and 
the type of each element, core or reflector. This has been 
simplified us the designation of all core elements by ITYPE=0 
and reflector elements by ITYPE=1. 

2. Programming Notes 

a. The property values in the two regions must be 
provided as an integral part of the program. 

b. the type of element, core or reflector, must be 
provided. 

c. The routine will process an unlimited number of 
models. It, therefore, must be halted when all data has been 
Merlized. This is done by specifying the final value of the 
numberwof elements in an IF-STOP statement. 

5. Parameters 

a. NEL - number of elements in the discretized model 

Det PE = type of element; ITYPE-0 is a core element 
and ITYPE=] is a reflector element. 

Gye Dees Vector Containing the diffusion length of the 
Gore and rezlectror elements. 

d. SGA - vector containing the absorption cross- 
section of the core and reflector elements. 

e. SGF - vector containing the fission cross-section 
Gf the core and reflector elements. 


£.. V= vector containing the neutron velocity. 
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g. ALPHA - vector containing the reactivity 
Pemperature coefficient. 


C. NODAL POINT COORDINATES AND ELEMENT 
NODAL POINT CONNECTIVITY 


i unpose 
| Nodal point positioning in the model is readily 

obtained in data deck form from this routine. Once the model 
to be investigated has been discretized into the finite ele- 
mene grid desired, the user provides the number of vertical 
and horizontal nodal points and their dimensional position 
along the axes. By the use of a nested loop, the nodal points 
are consecutively numbered and assigned the proper coordinate 
dimensions. 

ahevetemenemconnectivity, that is, the nodal points 
for each element, is also resolved by the use of a nested 
loop and several counters. Each iteration will yield two 
elements with their respective nodal point boundaries. This 
will continue until the nodal points for each element have 
been computed. 

iiewomteput Of thas rowtine will consist of two data 
decks. The first will provide the nodal point with its ver- 
tical and horizontal coordinates. The second will yield the 
element nodal point connectivity. 

Omitted from the element connectivity data deck, but 
required by the thesis program, is the type of element, core 


or reflector. This must be added to the deck individually. 


16 





Example of output (See Figure 1): 


nodal point coordinates 
nodal point horizontal position Vertreal position 


it 0.0 0.0 


element connectivity 
element nodal point connectivity type element* 


1 1 9 10 0 


*Note: must be added manually to data card. 


a Programming Notes 


a. The maximum number of vertical and horizontal 
nodal points in the discretized model must be provided. 

b. The horizontal and vertical positions of the 
discretization must be provided. 

c. The boundary points of the model will be numbered 
such that these nodal points are the last in the numerical 
Seauence. 

d. The routine is written such that it will process 
an indefinite number of models. Therefore, it must be halted 
when all the data has been utilized. This will be accomplished 
by specifying the maximum number of vertical points in the 
model in an IF-STOP statement. 

5. Parameters 

a. NH - maximum number of horizontal nodal points 
Mimole  dlsSeret1zed model . 

b. NV - maximum number of vertical nodal points in 


the discretized model. 
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c. X - vector containing the horizontal positions 
Gfethe NH points. 

d. Y - vector containing the vertical positions of 
Ene NV points. 

e. SYSNOD - vector providing the total number of 
nodal points. 

Pek “mVeCtor provadingmthe radial position of ithe 
nodal points. 

Go Spee vcctor providingethe gertical positdion of 
the nodal points. 

h. ELNOD - array providing the nodal point boundaries 
hoe each element, 

1. NEL - the total number of elements in the discre- 


tized model. 


D. NODAL NEIGHBOR CONNECTIVITY 
i Puspose 

This routine produces a data deck for the Crank- 
Nicolson and Implicit Gear methods containing each nodal 
point and the nodal points connected to it by the finite 
Emetiemn discretization. 

The regular rectangular grids are such that no more 
than six nodal points contribute to any one; therefore, an 
Nx7 array can be formulated for nodal neighbor connectivity. 
For example, nodal point 26 in the 112 element grid would be 
stored as follows (Figure 4): 


26 1 18 ay. 55 34 25 
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and nodal point 58 
58 49 50 59 S) 0 0. 

As written, zeros are inputted whenever a nodal point does 
not have six neighbors to maintain the symmetry of the matrix. 
This could, however, be compacted further by reading the data 
pmeaeoingle vector with a starting point counter to indicate 
when the next nodal point has been reached in the sequence. 
This would eliminate the need for the zeros to be supplied. 

This routine will provide a data deck for any 
rectangular grid and for an indefinite number of grids. 

2, weeloCoriehm 

This algorithm provides a dense, compacted matrix 
which reduces the storage size required by the main program. 
The use of a finite element method allows a certain amount of 
Compacting, i.e., the banded matrix. This is a function of 
the finite element size, that is, the Preeren the number of 
Meeal points, the larger the band. The size of the band is 
determined by the largest difference between nodal neighbors 
plus one. For example, in Figure 1, node ll has a maximum 
difference of 20 - 2 = 18, while at node 16 the difference 
is 43 - 7 = 36 which is the largest difference in the forty- 
five node system. As the number of nodal points increases, 
this maximum band width also increases. In Figure 2, the 
maximum difference is at node 16 (67 - 7 = 60); and, in 
Figure 3, the maximum difference of 124 - 10 = 114 exists 
at node 22. Therefore, although the size of the system is 


reduced from NxN to Nxq, the size is not constant and may 


ie 





not be small. In general, if one used banded storage, the 
numbering would proceed consecutively in the direction of 
fewest nodes. However, this would require the identification 
of those nodes on the boundary since they are of constant 
value and do not require integration. 

In the particular scheme of discretization utilized 
in this project, no nodal point has more than six nodal point 
contributors. This allows the matrix to be compacted further 
from Nxq to Nx/7. 

3. Programming Notes 

a. The total number of nodal points, the maximum 
number of nodal point contributors and the maximum number of 
horizontal and vertical nodal points must be provided. 

b. The routine will satisfy any rectangular grid. 
However, it must be instructed when a sufficient number of 
nodal points have been generated. This is accomplished by 
using the total number of nodal points desired in an [F-GO TO 
statement. 

c. This routine requires a positive means of stopping. 
This is accomplished by the use of the maximum number of 
vertical nodal points in the last data set in an IF-STOP 
Statement. 

GeameiiemOoutput Of this 7outine will place the central 
Modal point first with the contributors following. 

4. Parameters 
a. NV - number of vertical nodal points in the model. 


b. NH - number of horizontal nodal points in the model. 
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¢. NVH - total number of nodal points in the model. 
NVH = NV x NH. 

d. LCON - maximum number of contributing nodal points. 
bECON = 7). 

e. MNOD - an array, NVH x LCON, providing the central 


and contributing nodal points. 
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III. CRANK-NICOLSON METHOD OF SOLUTION 


Ae DESCRIPTION 
The Crank-Nicolson formulation is a numerical method 
Originally presented by J. Crank and P. Nicolson in 1947. 
Crank-Nicolson, being an implicit method, does not require 
the inverse of the matrix to be calculated. Therefore, 
advantage is taken of the sparse matrix inherently provided 
by the finite element method. If the inverse of the sparse 
Meerix 1s formed, it will be full. 
The Crank-Nicolson method will solve the set of linear 
differential equations represented by 
N 
hs 


N 
a ew) elena e 
jh Ge ages aN (3) 


A.. w.(t) = 
1) v5 Ct) ieee’ 


Crank-Nicolson, after some algebra, yields 


Bee) | ee EAT (4) 
At Z 

The implicit system (4) may be solved by an iterative process; 
in this case, the Gauss-Seidel method is used to solve the 
set of simultaneous algebraic equations (4) formulated. Of 
all the stable numerical methods in the case of single step 
implicit, Crank-Nicolson has the smallest truncation error [2]. 
The Gauss-Seidel method of iteration is continuously using the 
newest solution values allowing convergence to the solution 


to be the most rapid. 
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BeeeeebeCl OF DEGREES OF FREEDOM 

As might be expected, as the grid size became finer, the 
computer core requirements increase. In addition, more time 
is required for the problem to reach a steady state solution. 
At the same time, as shown in Figures 5-7, the solution values 
converged, as the mesh size became finer, to a more accurate 
solution. 

As shown in Table I and Figure 8, storage and processing 
time, that is, CPU time to reach a steady state solution, in- 
creased markedly between a forty-five node grid and a one 
hundred thirty-two node grid. This yielded a fifty-nine 
PeGeent increase 1n storage and better than a five hundred 
Pescemc increase In processing tame. At the same time there 
is an improvement in the solution of only four and one-half 
pemcent . 

Comparison to a seventy-two node grid yielded far more 
Satisfying results. At seventy-two nodes, there is an eighty 
percent inerease in CPU time and a sixteen percent increase 
in storage requirements while obtaining a two and one-half 
Demecent Amprovement in the solution. 

Similar results were obtained in problem time to solution. 
For one hundred thirty-two nodes, the time more than doubles; 
while for seventy-two nodes, the time was increased by twenty 


percent. 


fee rrect OF ERROR CRITERION 
The error criterion used throughout this thesis informs 


the integration routine when it has achieved sufficient 
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Seceement between iterates (i.e., when it can stop and go 
Bomthe Next iteration step). 

In Crank-Nicolson, the error:criterion simply consisted 
of the difference in successive Gauss-Seidel iterates divided 


by the current iterative value. 


it 4 


jie error criterion was varied from 10 ~ to 10 ° to 
observe the effects on the solution and computer requirements. 
For Crank-Nicolson the steady state solution value remained 
the same regardless of the error criterion. This demonstrated 
miawmmne ECXtra time required to satisfy a tighter error cri- 
terion at each ee anion was not necessary to reach a satis- 
factory steady state solution. As was expected, the computer 
processing time did increase significantly as Eieeunare 
criterion was made smaller. An unusual result occurred when 
Bae E€xror Criterion was set to oe. This particular value 
resulted in a processing time less than. that required for 
‘ge It would appear that this might have been the result 

Ge two contributing factors. With the closeness of the error, 
each time an iteration was performed, it was using a better 
solution, and each new time step also had a better solution 
value. Although there were many more time steps required 

for this error criterion, the steady state value was rapidly 
approached because fewer iterations were required at each new 
fee.  tnese results are shown in Table II. These results 


: do not imply 


peineasericciterion between 10 -.and 10 
that this would be the result obtained for all problems of 


this type. 
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Ueeeetrecl OF INCREASE IN TIME STEP SIZE 

In this method, the time step can be increased or not 
depending on the solution. If a solution is not obtained 
with a particular time step, that same time step is reduced; 
and the routine attempts to attain solution with the smaller 
time step. This continues until either a satisfactory solu- 
tion value is obtained or the built in default value is 
reached which terminates the routine. When a solution is 
attained, the routine compares the number of iterations 
required to reach that value to that of the previous time 
step and to a specified number of iterations. If the current 
number of iterations is less than either of those numbers, 
the program increases the time step by an input constant 
value. 

The change made to the initial time step as the steady 
state was approached was the critical value if the solution 
was to converge to a steady state. It appeared that Crank- 
Nicolson was particularly sensitive to the amount the time 
step was increased as the solution was approached (Figure 9). 
The method was unable to recover if, when nearing the knee of 
the curve, the size of the increase was such that the steady 
state value was exceeded. Once the solution was passed, the 
results indicated a diverging oscillation about the steady 
state value. 

Several values of the initial time step were utilized in 
all grids to attempt to locate the steady state. The results 


Biethese trials are provided in Table III. {It is noteworthy 


Zo 





that if the increase was greater than 1.2 times the initial 
step size, the steady state would never be achieved. Whether 
Piss would be the case for all grids is not obvious. This is 
particularly evident in view of the results obtained by 

Olsen [6], who utilized an increase of 1.5 times the initial 
Seep tor a thirty-eight node grid. It would, therefore, seem 
that a trial and error approach would be necessary until an 


increase in step size resulted in a steady state solution. 
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IV. DVOGER (GEAR) METHOD OF SOLUTION 


ev DESCRIPTION 

The DVOGER (Gear) routine is an IMSL library routine 
which integrates a system of explicit first order differ- 
ential equations. Within this library routine is the capa- 
bility of solving both stiff and nonstiff systems by selection 
of an indicator. In solving the nonstiff system, the Adams 
predictor-corrector is utilized. For the stiff system, Gear's 
predictor-corrector is used. Gear's method computes the 
Jacobian for the system of ordinary differential equations 


Mimonder tO Optimize the time step for each integration. 


Bee EPFECY OF DEGREES OF FREEDOM 

Gear's method placed a much greater demand on the computer 
in terms of processing time and storage than the other methods 
considered. This was primarily due to three causes: 1) the 
calculation of the Jacobian, 2) the transformation of Equation 
Wo tovexplicit form, and 3) the initial absolute value of 
the solution. The first two items require a substantial 
increase in both CPU time and computer storage. Several 
variances were implemented in this program to determine if 
better use could be made of the method to make it competitive 
with Crank-Nicolson or the Implicit Gear methods. 

When the system was considered nonstiff, the storage 
requirements decreased since the Jacobian was not calculated; 


but, the processing time increased by one order of magnitude 


oy 





when compared to the stiff system treatment due to the small 
Em@meesceps' taken to solution. The routine calls for an 
initial absolute solution value of one to be supplied. This, 
however, adversely affects the progress of the problem since 
this value of one is initially compared to values of Tas 
iis limits the time step taken by the routine initially and 
continues to affect the progress until the steady state solu- 


tion is neared. When an initial value of 1014 


was utilized, 
the processing time was decreased by forty percent. This was 
due to the fact that larger time steps were allowed from the 
Oieset, Since the previous solution value (in this case the 
initial value) was comparable to the value obtained in the 
first calculation. 

At best, when utilizing Gear's method, the processing 
times were two orders of magnitude greater than that required 
for the Implicit Gear method and twenty times the Crank- 
Nicolson processing times for the 132 DOF system, as shown in 
Table I. Core requirements were also increased significantly, 
particularly at one hundred thirty-two degrees of freedom. 
For this grid, DVOGER (Gear) was approximately three times 
Crank-Nicolson and the Implicit Gear core requirements. 

Gear's method provided the same rea crate solution 
values as the other two methods, and the transient curves 


were very similar to those obtained in the Implicit Gear 


method (Figures 10-15). 
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SemeeerreECc! OF ERROR CRITERION 

ineGear s method, the time step size is adjusted so that 
the single step error estimate divided by the previous maximum 
sOlution value is less than the error criterion in the 
Euclidean norm. The single step error estimate is a multiple 
of the difference between the predicted and corrected values 
of the variable. 

The error criterion, Table II, was varied for this method 
as it was for the others. The results were the same; as the 
criterion was eeenecr CPU time increased. In this case, 
however, the time became excessive very rapidly, more than an 


eee ge eee ao oeeeehicuirs far) 10s for che model 
with 45 DOF. Similarly, the solution values were the same 


regardless of the criterion utilized. 
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V. IMPLICIT GEAR METHOD OF SOLUTION 


Pa DESCRIPTION 

This numerical method is particularly useful for the 
solution of large, sparse systems of implicit, stiff differ- 
ential equations [4]. In contrast to Gear's method, the 
Implicit Gear method eliminates the need of an exact NxN 
Jacobian, but rather requires only an exact Nx7 Jacobian. 
This is due to the fact that Equations (3) are handled directly, 
thereby retaining the sparseness of the matrix. 

Gear's predictor-corrector method and the compact matrix 
makes use of storage because the implicit Equations (3) are 
handled directly. This results in very efficient use of the 
computer both in terms of storage and processing time require- 
ments. 

The user is required to provide a subroutine that evalu- 
ates the system of equations being investigated, as well as, 
for efficiency, a subroutine to evaluate the Jacobian of the 
system of equations. 

This program is a modification by Franke [4] to the 
routine DFASUB [7]. One of the major changes to the routine 
ice riemerecatment Of the error. In the Implicit Gear 
Method, instead of using the Euclidean norm, the root mean 
Square norm is utilized. This is no more than the Euclidean 
norm divided by the square root of the number of components. 
In addition, the maximum value of the component is updated 


before the norm of the relative error is computed. 
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Pee EPPECT OF DEGREES OF FREEDOM 

The steady state values obtained in this method were the 
same as those obtained in Crank-Nicolson and DVOGER (Gear) 
as shown in Figures 11-18. This result does not imply that 
these methods will always give identical steady state solu- 
tions. Certainly they do give different transient solutions. 
The major advantages of this method are 1) the computer 
processing time was reduced as the grid size became finer 
ome Z)jethe storage requirements are slightly greater than 
that required in the Crank-Nicolson method due to the compu- 
tation of the Jacobian and the storage of up to seven past 
solutions which control the size of the time step. Implicit 
Gear requires about NxZ5 more storage locations than Crank- 
Nicolson. This amounts to about only thirteen thousand bytes 
for the one hundred thirty-two degrees of freedom system 
(single precision). The magnitudes of the increase in pro- 
cessing time dropped dramatically in this method. For the 
Same initial disturbance, the time increased eighty percent 
for the one hundred thirty-two node grid and twenty-four 
percent for the seventy-two degrees of freedom. In addition, 
Implicit Gear gave a more accurate transient solution. The 
results of this method are shown in Table I and comparison 


to the other methods will be made in Chapter VI. 


fee EEEECI OF ERROR CRITERION 
MimimperettrGear, the error criterion is defined as in 


DVOGER (Gear) except that the root mean square of the Euclidean 


Sl 


* 
ee 





norm and the updated maximum solution value are used in the 
computation. 

The error criterion was varied substantially to determine 
what effect it had on computer storage and time requirements, 
feewelleas its effect on the accurdcy of solution. As shown 
in Table II and Figures 19-20, a wide range of convergence, 
iO to Hon, was investigated with interesting results. 
Expectedly, the processing time did increase as a tighter 
error criterion was required. However, the criterion had 
little effect on the track through the transient solution for 
values of less than a When using an error of one-half 
and one, the transient solution varied substantially. The 
same steady state solution was obtained, although at a later 
point in problem time as shown in Figure 19. For the stiffer 
error criterion requirements, the processing time necessary 
is doubled by utilizing a criterion of 10°* instead of a value 
it 


greater than 10 This, certainly, did not appear to warrant 


the closer tolerance level. 
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Vee REOULTS 


A. COMPARISON OF TIME AND STORAGE REQUIREMENTS 

With all three methods yielding the same steady state 
solution (Figures 10-18), the comparison of methods became 
one of time and storage requirements placed on the computer 
rather than one based on which provided the best solution. 
In the transient stage, Gear tracked slightly better than 
Implicit Gear, but both Gear and Implicit Gear tracked better 
than Crank-Nicolson (See Figures 10-18). For the three systems 
of equations utilized, the Implicit Gear method performed 
significantly better than either Crank-Nicolson or DVOGER 
(Gear) in CPU time and was comparable to Crank-Nicolson and 
superior to Gear in storage requirements (Table I). This 
becomes more significant as the number of degrees of freedom 
increase. The Implicit Gear is less sensitive to the increase 
in size of the system of equations as shown in Figure 8. 

Comparing the forty-five and the one hundred thirty-two 
degrees of freedom, the core requirements increased slightly 
and the time doubled for the Implicit Gear method. In con- 
trast, for Crank-Nicolson, the time increased by six times 
and the core by 100K bytes; and, for DVOGER (Gear), there 
was a ten times and 500K byte increase in time and storage 


hespectively . 
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Bee DEGREE OF FREEDOM EFFECT ON SOLUTION 

Varying the degrees of freedom resulted in a better 
solution value as shown in Figures 5-7 and 21-26. There was 
a four and one-half percent better solution obtained for the 
one hundred thirty-two degrees of freedom compared to the 
forty-five degrees of freedom. This better solution is very 
costly both in terms of computer core requirements and CPU 


time (Figure 8). 


c ERROR CRITERION 

Error criterion was varied for all methods to determine 
the effect on solution and computer. As might be expected, 
the computer processing time increased for all methods as 
more rigid tolerance was imposed (Table II). However, although 
the criterion was made very close, there was no appreciable 
effect on the transient solution until the error criterion 
was greater than fae and no effect on the steady state 
Solution. This is significant in that, as stated above, the 
Preecessing time increases greatly as the criterion is tight- 
ened. This shows that for the particular problem considered 
here, some close error criteria yield no advantage, only the 


disadvantage of requiring more time in the computer. 


D- INITIAL DISTURBANCES 

Tiewenpat Of Various initial disturbances (central, skewed 
at the core-reflector interface and uniform throughout the 
core) yielded the same steady state solution value. The 


transient solution varied substantially due to the input 
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position of the disturbance. In most cases, due to the 
Magnitude and scope of the initial disturbance, the uniform 
disturbance required the least processing time to reach the 


steady state value. 


E. COMPARISON OF SOLUTION TRACKING 

Two nodal points were investigated graphically for each 
disturbance, integration method and DOF, as shown in Figures 
5-7, 10-18 and 21-26. At these test points, the methods 
performed similarly in all cases in both the transient and 


steady state solutions. 


oe) 


Vit. SiEST PROBLEMS 


fee PHYSICAL MODEL 

A cylindrical reactor with the dimensions and properties 
given in Table IV was used as a model. A radial slice of 
this model was discretized by various finite element grids 


as shown in Figures 1-3. 


B. COMPUTER PROCESSING CONSIDERATIONS 

The programs presented in this thesis, Appendix A, were 
written in the FORTRAN IV language and all computer runs, 
Table V, processed on the IBM 360/67 computer using the 
FORTRAN 'G' compiler. It is expected that the FORTRAN 'H' 
compiler would have supplied results similar to those obtained 
with the FORTRAN 'G' compiler. This was not done because the 
FORTRAN 'H' compiler requires 350K bytes as a minimum core 
requirement. The majority of the runs performed in this 
thesis required significantly less than 300K storage. Single 
precision (six to seven significant digits) was used through- 
out this research. As has been shown [6], double precision 
solutions are at variance by less than 0.01 percent with 


single precision solutions. 


G> PROBLEM ANALYSIS 
The techniques used by Olsen [6], modified as described, 
and the Implicit Gear method [4] were utilized to solve the 


problem delineated in Ref. [5]. The problem was approached 
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using finite element discretizations with forty-five, seventy- 
two and one hundred thirty-two degrees of freedom and providing 
an initial disturbance. Three initial disturbances were pro- 
vided as follows: central at the core center (R = 0 cn., 

Z = 0 cm.); uniform throughout the core; and, a skewed dis- 
Eumbamce at the core-reflector interface (R = 60 cm., Z = 0 


em.). 


D. PROGRAM USAGE 

In order to properly utilize the routines presented in 
this thesis, several points must be considered. 

In the Crank-Nicolson routine, there was difficulty in 
obtaining a steady state solution. The ultimate cause was 
the size with which the previous time step increases. This 
requires a trial and error approach for the particular grid 
eze. For this project, an increase of ten or twenty percent 
feteds Sdtistfactory results; but, anything larger results in 
a divergent oscillation about the steady state value. 

In all three methods, the dimensioning should be set by 
the particular problem being solved, i.e., determined by the 
degrees of freedom. If this is done, most efficient use will 
be made of the computer, and the user will experience better 
turn around time on the equipment. 

In the use of the Implicit Gear method, the data cards 
for nodal neighbor connectivity must have the contributing 
nodes listed consecutively, with the zeros, used for nodes 
not having six neighbors, being last on the cards. For example, 


in the one hundred thirty-two degrees of freedom system 


oF 





(Figure 3), the nodal neighbor connectivity for node 122 


would be written as 


iZ2 oS et 0 0 0 We 


The data deck arrangement for all three methods is given in 


Table VI. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

It has been shown in this research project that of the 
three methods investigated in the solution of the test problem 
that the Implicit Gear performed in a superior manner in CPU 
time requirements and was comparable to the storage require- 
ments of Crank-Nicolson. In addition, Implicit Gear was the 
least sensitive to the change in the degrees of freedom. 

The method used had no effect on the steady state solution 
value obtained as all three yielded the same results for the 
same DOF. In the transient solution, Implicit Gear conformed 
very closely to DVOGER (Gear). As shown in Figures 10-15, as 
the DOF increased, the transient solutions moved in the direc- 
tion of the transient solution of Gear's method. 

As the DOF in the discretized finite element was increased, 
a better solution was obtained. This, however, was counter- 
fered by the fact that for the small increase in accuracy 
obtained by the finer mesh, the time and storage need of the 
computer increased significantly. 

A very important result of this project was the peices of 
varying the error requirements of the problem. This yielded 
essentially no change in the accuracy of the transient solu- 


i 4 


to 10° and no variance in the steady 


4 


tion for a range of 10. 
Seaee solution for a range from 1.0 to 10 In using these 


methods of integration, the error criterion selected would 
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eemoecunction of the solution desired. If the user is only 
concerned with the steady state solution, the tolerance could 
be relaxed to 1.0. However, if the transient solution is 
needed, an error of Om appears to be the least rigid value 


which will still provide a satisfactory track to steady state. 


B. RECOMMENDATIONS 

It would be of value to investigate further the Implicit 
Gear method through the use of other problems solving a non- 
linear system of ordinary differential equations. In addition, 
larger degrees of freedom should be attempted to evaluate any 


limitations that may exist in this method of integration. 
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TABLE I 


COMPARISON OF SOLUTION METHODS 


METHOD 


CRANKO 
DVOGER (GEAR) 


IMPLICIT GEAR 


CRANKO 
DVOGER (GEAR) 
IMPLICIT GEAR 


CRANKO 
DVOGER (GEAR) 
IMPLICIT GEAR 


Processing time is that required to reach a steady state 


solution. 


CORE 


158 
148 
124 


184 
244 
124 


610 
124 


(K) 


Central 


PROCESSING TIME* (min) 


45 Nodes 


Od 
Lye 


05 


72 Nodes 


89.5 


132 Nodes 


41 


9.1 
>400 
1.9 


Skewed 


oe 


114.1 


9.1 
>400 


Uniform 


t. 25 


Oru 
>400 





TABLE II 


EFFECTS OF SOLUTION ERROR CRITERION 


45 Nodes 
SOLUTION ERROR PROCESSING TIME TIME AT SOLUTION 
CRITERION (min. ) (sec. ) 


CRANK-NICOLSON 


0.1 1.8 Soon me 
0.01 55 PTO Oe 
0.001 14.2 ese lon. 
0.0001 0.83 Beene 
IMPLICIT GEAR 
1.0 0.99 7S ones 0 ae 
nes 0.92 ig alan @ 
0.1 1.05 7aee 0 = 
0.01 1.20 7.64 x 107° 
0.001 1.45 sieio x 1m 
0.0001 1.58 6.44 x 107° 
DVOGER (GEAR) 
0.1 18.6 5.64 x 10> 
0.01 29.5 5.79 x 10° 
0.001 78.8 es ene 
0.0001 >300 = 
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TABLE III 


Peeeers OF VARYING THE INCREASE OF THE INITIAL TIME STEP 
IN CRANK-NICOLSON METHOD 


INCREASE CORE PROCESSING TIME SOLUTION TIME 
(min. ) (sec.) 
45 NODES 
1), 158K 4.99 acme 
1.2 158K Le arate 
aes 184K 7.42 OSCILLATING 
72 NODES 
eal 252K as 4.38 x 10°° 
Le 184K 3.05 3.23 x 10> 
iS 252K as OSCILLATING 
1.4 252K 9.1 OSCILLATING 
no 184K eel OSCILLATING 
132 NODES 
ae 
ped 252K onal 5.58 x 10 
lag 252K 8.4 8.03 x 10° 
ies 184K me OSCILLATING 
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SYMBOL 


coc kt WD DW 


< 


ac 


ar 


COMPUTER 
SYMBOL 


coi 2 7 1G 

R 

R 

Zz 

Zz 
V/VELOCT 
D/DSUBF 


D/DSUBC 


SGA/SGMAAF 


SGA/SGMAAC 


ZNU 


SGF/SGMAF 


SGE/ - 


FISFAC 
HBAR 


ALPHA 


TABLE IV 


PHYSICAL CONSTANTS 


DEFINITION 


Total radius 

Core radius 

Gore nerlene 

ioral heroic 
Neutron velocity 
Neutron diffusion 


coefficient (core) 


Neutron diffusion 
coefficient 
(reflector) 


Neutron absorption 
GEoss,- Section 
(core) 


Neutron absorption 
cross-section 
(reflector) 


Number of neutrons 
per fission 


Neutron fission 
cross-section 
Ceome)) 


Neutron fission 
cross-section 


Fission ehergy 


Modified convection 0.0632 cal/em” 


heat transfer 
coefficient 


Reactivity temper- 
ature coefficient 
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VALUE 


90 cm. 


Og) em, 


PoUeen. 


220 ene 


4 
0 


Te 


it 


“OO smem | 


aeetae 
~915 Cm. 


cine Ic 


.01401 cm 


652x10 


ine 


1 


1 


tz 


/ °C 


cm/sec 


Caly/-ers 


-sec-°C 





TABLE V 


LIST OF COMPUTER RUNS 


LS 


RUN METHOD NODES DISTURBANCE COVERGENCE TIME STEP 
CHANGE 
1 CRANK-NICOLSON 45 CENTRAL Ol Pei 
2 ieee 
3 1.5 
4 0.01 eye 
5 0.001 
6 0.0001 
7 UNIFORM Oral 
8 SKEWED 
9 qe CENTRAL 0.1 10 
10 a 
et rez 
12 ies 
13 iA 
14 1.5 
15 UNIFORM eee 
16 SKEWED 
17 132 CENTRAL 0.1 dean 
i? 
19 1.5 
20 UNIFORM 122 
21 SKEWED 
22 DVOGER 45 CENTRAL 0.1 N.A. 
MTH=0 
23 DVOGER Det 
MTH=1 
YMAX=1 
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TABLE V (Continued) 


RUN 


24 
25 
26 


27 
28 
29 
30 
31 
BZ 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 


47 


METHOD NODES 
DVOGER 45 
MTH=1 
YMAX=1 
V2 
Loi2 


IMPLICIT GEAR 45 


Le 


IES Z 


45 


CENTRAL 


UNIFORM 
SKEWED 


CENTRAL 


UNIFORM 
SKEWED 
CENTRAL 
UNIFORM 
SKEWED 
CENTRAL 


UNIFORM 
SKEWED 
CENTRAL 
UNIFORM 
SKEWED 
CENTRAL 
UNIFORM 
SKEWED 
CENTRAL 
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DISTURBANCE CONVERGENCE TIME STEP 


CHANGE 





faeLeE V (Continued) 


RUN METHOD NODES DISTURBANCE CONVERGENCE TIME STEP 
CHANGE 
48 a 
49 DVOGER 45 CENTRAL Os 
MTH=1 14 
YMAX=10 
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TABLE VI 
DATA DECK ARRANGEMENT 


CRANK-NICOLSON 
prt le 
NUMEL ,NUPBP ,NUMSNP, NFULEL 
List of outer boundary points 
MTH , MAXDER,NCOUNT 
ZNU,FISFAC,HBAR, EPSVAL, ERRVAL ,AFUEL 
TO,H,TF,HMIN ,HMAX 
Vo oeieutron velocity 
D - diffusion length 
oGA - absorption cross-section 
SGF - fission cross-section 
ALPHA - reactivity temperature coefficient 
PSIIV - initial disturbance flux 
System nodal points with axial and radial position 
Element nodal connectivity 


Nodal neighbor connectivity 


DVOGER 
Title 
NUMEL , NUPBP , NUMSNP,NFULEL 
List of outer boundary points 
MTH ,MAXDER, NCOUNT 
ZNU,FISFAC ,HBAR, EPSVAL , ERRVAL, AFUEL 
TO,H,TF,HMIN , HMAX 
V 
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TABLE VI (Continued) 


D 

SGA 

SGF 

ALPHA 

PSIIV 

System nodal points with radial and axial position 


Element nodal connectivity 


IMPLICIT GEAR 

NCORE 
with uniform initial disturbance 

List of core elements 
Tetle 
NUMEL,NBP,NUMSNP , NFULEL,NELROW ,MXEVAL, NCOUNT , NCMPK 
NDE,NL 
List of outer boundary points 
VELOCT , DSUBF, DSUBC, SGMAAF , SGMAAC , SGMAF ,, FISFAC ,HBAR 
ZNU , ALPHA ,RMSEPS,AFUEL,TSTART,TEND,PINITV,SCALE 
HMIN , HMAX 
MTH ,MAXDER,NPROB,NTYPE,JPLOT ,NPLOT, IRUN 
Nodal neighbor connectivity 


System nodal points with radial and axial position 


Element nodal connectivity 
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Figure 4. Nodal Neighbor Connectivity 
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C-----------NUCLEAR PRIPERTY DATA GENERATOR------------ 
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OR TO PROVIDE NODAL NEIGHBOR 
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